{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f414c9c0-89d1-4408-a658-68de4c05bb25",
   "metadata": {},
   "source": [
    "# Cholesky decomposition\n",
    "\n",
    "A [Cholesky decomposition](https://en.wikipedia.org/wiki/Cholesky_decomposition) is a useful factorization of Hermitian, positive-definite matrices into the product of a lower triangular matrix $L$ with its conjugate transpose $L^{*}$.\n",
    "\n",
    "Numpy has a function [numpy.linalg.cholesky](https://numpy.org/doc/stable/reference/generated/numpy.linalg.cholesky.html) built-in for computing Cholesky decompositions. Cunumeric also implements this function, and it can be used as an immediate drop-in replacement.\n",
    "\n",
    "<details>\n",
    "<summary>License</summary>\n",
    "<pre>\n",
    "\n",
    "Copyright 2023 NVIDIA Corporation\n",
    "\n",
    "Licensed under the Apache License, Version 2.0 (the \"License\");\n",
    "you may not use this file except in compliance with the License.\n",
    "You may obtain a copy of the License at\n",
    "\n",
    "     http://www.apache.org/licenses/LICENSE-2.0\n",
    "\n",
    "Unless required by applicable law or agreed to in writing, software\n",
    "distributed under the License is distributed on an \"AS IS\" BASIS,\n",
    "WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n",
    "See the License for the specific language governing permissions and\n",
    "limitations under the License.\n",
    "</pre>\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "389cd191-ccda-4597-8e08-8d01ac226bee",
   "metadata": {},
   "source": [
    "To get started, `import cunumeric as np` (just the same way we would import `numpy`)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8bb45fad-6f50-40de-ac68-083aa5fe0f1b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import cunumeric as np  # instead of numpy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ef2bc57-e703-40ce-8aaa-d45408259c7a",
   "metadata": {},
   "source": [
    "At this point we can call `np.linalg.cholesky`, exactly how we would with Numpy, but will get the result computed by Cunumeric's `cholesky` function. Let's quickly try it out with a simple identitity matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "2933b8cb-406c-4ebd-999e-4f3f3c687a6f",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 0., 0., ..., 0., 0., 0.],\n",
       "       [0., 1., 0., ..., 0., 0., 0.],\n",
       "       [0., 0., 1., ..., 0., 0., 0.],\n",
       "       ...,\n",
       "       [0., 0., 0., ..., 1., 0., 0.],\n",
       "       [0., 0., 0., ..., 0., 1., 0.],\n",
       "       [0., 0., 0., ..., 0., 0., 1.]])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.cholesky(np.eye(100))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09b58cf2-b7ee-45eb-ade6-4d99d930047b",
   "metadata": {
    "tags": []
   },
   "source": [
    "We'd like to get some information about how well Cunumeric's `cholesky` function performs. In order to obtain accurate timings, we need to use the `time` function from `legate.timing`. Let's define a helper function `cholesky_timed` that calls the `time` function for us, and prints out the results as well:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "a5af9db1-3656-447f-887e-e08dc0f287c8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Because of Legate's deferred execution model, legate.timing should be used instead\n",
    "# of standard Python datetime utilities. Python datetime.now would return the time\n",
    "# a task is *scheduled*, not necessarily the time a task finishes executing.\n",
    "from legate.timing import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ea1027df-b4b5-46d8-a1a0-73ba8c0a1c73",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def cholesky_timed(A):\n",
    "    start = time()\n",
    "    result = np.linalg.cholesky(A)\n",
    "    stop = time()\n",
    "    \n",
    "    n = A.shape[0]\n",
    "    flops = (n**3) / 3 + 2 * n / 3\n",
    "    total = (stop - start) / 1000.0\n",
    "    print(f\"Elapsed Time: {total} ms ({(flops / total):0.2} GOP/s)\")\n",
    "    \n",
    "    return result"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ebbfeba-c768-47b9-b088-32de999b051f",
   "metadata": {},
   "source": [
    "Now we can define a matrix $A$ to decompose. Let's make some thing a little more interesting than a plain identity matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0678173d-1eb6-4ad6-b282-9b41e9cb3b40",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "A = np.eye(1000, dtype=np.float64)\n",
    "A[:, 2] = 1\n",
    "A = np.dot(A, A.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "06a767f0-b7d7-4eec-8ad8-8d1c127babff",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Elapsed Time: 31.277 ms (1.1e+07 GOP/s)\n"
     ]
    }
   ],
   "source": [
    "L = cholesky_timed(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6507a3d8-9667-4fac-bf50-f36241f4b11b",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
